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We make a systematic study of the cosmological dynamics for a number of f(R) gravity theories 
in Palatini formalism, using phase space analysis as well as numerical simulations. Considering 
homogeneous and isotropic models, we find a number of interesting results: (i) models based on 
theories of the type (a) f(R) = R - /3/R n and (b) f(R) = R + am R - j3, unlike the metric 
formalism, are capable of producing the sequence of radiation-dominated, matter-dominated and 
de-Sitter periods, and (ii) models based on theories of the type (c) f(R) = R + aR m — /3/R n can 
produce early as well as late accelerating phases but an early inflationary epoch does not seem to 
be compatible with the presence of a subsequent radiation-dominated era. Thus for the classes of 
models considered here, we have been unable to find the sequence of all four dynamical epochs 
required to account for the complete cosmological dynamics, even though three out of four phases 
are possible. 

We also place observational constraints on these models using the recently released supernovae 
data by the Supernova Legacy Survey as well as the baryon acoustic oscillation peak in the SDSS 
luminous red galaxy sample and the CMB shift parameter. The best-fit values are found to be 
n = 0.027, a — 4.63 for the models based on (a) and a = 0.11, j3 = 4.62 for the models based on (b), 
neither of which are significantly preferred over the ACDM model. Moreover, the logarithmic term 
alone is not capable of explaining the late acceleration. The models based on (c) are also consistent 
with the data with suitable choices of their parameters. 

We also find that some of the models for which the radiation-dominated epoch is absent prior to 
the matter-dominated era also fit the data. The reason for this apparent contradiction is that the 
combination of the data considered here does not place stringent enough constraints on the cosmo- 
logical evolution prior to the decoupling epoch, which highlights the importance of our combined 
theoretical-observational approach to constrain models. 



I. INTRODUCTION 

Recent high-precision observations by the Wilkinson Microwave Anisotropy Probe (WMAP), together with other 
Cosmic Microwave Background (CMB) and high redshift surveys have produced a wealth of information regarding 
the early universe. The analysis of the resulting data has provided strong evidence for the core predictions of the 
inflationary cosmology, including the almost spatial flatness of the universe [J 0. Furthermore, these observations 
coupled with the low redshift supernovae surveys [!, 0, IE IS and observations of large scale structure d, Q and 
baryon acoustic oscillations [13] suggest that the universe is at present undergoing a phase of accelerated expansion (see 
Refs. [ll],[l2j for reviews). Thus a 'standard' model of cosmology has emerged which is characterised by four distinct 
phases: accelerated expansions at both early and late times, mediated by radiation-dominated and matter-dominated 
eras. The central question in cosmology at present is, therefore, how to successfully account for these distinct 
dynamical phases in the history of the universe, and in particular whether they can all be realised simultaneously 
within a single theoretical framework, motivated by a fundamental theory of quantum gravity. 

A number of scenarios have been proposed to account for these dynamical modes of behaviour. These fall into 
two categories: (i) those involving the introduction of exotic matter sources, and (ii) those introducing changes to 
the gravitational sector of General Relativity. Among the latter are f(R) gravity theories, which involve non-linear 
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generalisations to the (linear) Hilbert action. Nonlinear modifications are expected to be present in the effective action 
of the gravitational field when string/M-theory corrections are considered QS Q OS 0]| • 

A great deal of effort has recently gone into the study of such theories. An important reason for this interest has 
been the demonstration that generalised Lagrangians of this type - which include negative as well as positive powers 
of the curvature scalar - can lead to accelerating phases both at early [17| and late [l8|, |l9| times in the history of the 
universe (see also Ref. |20||). 

In deriving the Einstein field equations from the Hilbert action the variations are taken with respect to the metric 
coefficients, while the connections are assumed to be the Christoffel symbols defined in terms of the metric. An 
alternative procedure - the so called Palatini approach, originally considered by Einstein himself - is to treat both the 
metric and the (affine) connections as independent variables and perform the variations with respect to both. In the 
case of linear Hilbert action both approaches produce identical results, as long as the energy-momentum tensor does 
not depend on the connection [2l[. This, however, is not the case once the gravitational Lagrangian is allowed to be 
nonlinear. In that case the two methods of variation produce different field equations with nontrivial differences in 
the resulting dynamics. 

Performing variations of nonlinear actions using the metric approach results in field equations that are fourth order; 
which makes them difficult to deal with in practice. Furthermore, within this framework, models based on theories 
of the typ e f(R) = R — /3/R n have difficulties passing the solar system tests [22[ and having the correct Newtonian 
limit [23]. In addition such theories suffer gravitational instabilities as discussed in Ref. [2J]. Also recent studies 
have found that these theories are not able to produce a standard matter-dominated era followed by an accelerated 
expansion (2f| [26j. Finally, models based on theories of the type f(R) = R + aR' n — f3/R n have been shown to have 
difficulties in satisfying the set of constraints coming from early and late-time acceleration, Big Bang Nucleosynthesis 
and fifth- force experiments (27[. (See Refs. 2^1 for other works concerning the metric formalism.) 

Variations using the Palatini approach [29ll30l|. however, result in second order field equations which are in addition 
known to be free of instabilities [3ll 32| . Such theories have also been shown to satisfy the solar system tests and 
have the correct Newtonian limit [33 1. Here we shall concentrate on the Palatini approach and consider a number 
of families of f(R) theories recently put forward in the literature. These theories have been the focus of great deal 
of interest recently with a number of studies attempting to determine their viability as cosmological theories, both 
theoretically [HI, yH, [3f| [H, [37j and observationally [3a |. Despite these efforts, it is still fair to say that cosmological 
dynamics of models based on such theories is not fully understood. Of particular interest is the number of dynamical 
phases such models are capable of admitting, among the four phases required for strict cosmological viability, namely 
early and late accelerating phases mediated by radiation-dominated and matter-dominated epochs. Especially it is 
important to determine whether they are capable of allowing all four phases required for cosmological evolution. We 
should note that the idea that such theories should be capable of successfully accounting for all four phases is clearly 
a maximal demand. It would still be of great interest if such theories could successfully account for a sequence of 
such phases, such as the first or the last three phases. For example the presence of the last 3 phases could distinguish 
these from the corresponding theories based on the metric formalism j25[ . Also observational constraints have so far 
only been obtained for the models based on the theories of the type f(R) = R — (3/R n . Such constraints need to be 
also studied for other theories considered in this context in the literature. 

Here in order to determine the cosmological viability of models based on f(R) theories in Palatini formalism we shall 
employ a two pronged approach. Considering homogeneous and isotropic settings, we shall first provide autonomous 
equations applicable to any f(R) gravity theory. We shall then make a systematic and detailed phase space analysis of 
a number of families of f(R) theories. This provides a clear understanding of the dynamical modes of behaviour which 
are admitted by these theories, and in particular whether these theories possess early and late accelerating phases 
which are mediated by the radiation-dominated and matter-dominated eras respectively. The existence of these phases 
is a necessary but not sufficient condition for cosmological viability of such theories. To be cosmologically viable, it 
is also necessary that the subset of parameters in these theories that allows these phases to exist are also compatible 
with observations 0|. We constrain these parameters for a number of families of f(R) theories using the data from 
recent observations, including recently released supernovae data by the Supernova Legacy Survey [7[ as well as the 
baryon acoustic oscillation peak in the Sloan Digital Sky Survey (SDSS) luminous red galaxy sample [To| and the 
CMB shift parameter 0,0). 

The plan of the paper is as follows. In Section [IT] we give a brief account of the Palatini formalism and provide the 
basic equations for general f(R) theories. We also introduce new variables which allow these equations to be written 
as autonomous dynamical systems. In Section IIIII we proceed to study the cosmological dynamics for a number of 
classes of f(R) theories. Particular emphasis will be placed on finding cases which can allow the largest number of 
dynamical phases required in the cosmological evolution. In Section IIVI we obtain observational constraints on these 
theories. This allows further constraints to be placed on the parameters of these theories. Finally we conclude in 
Section El 
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II. f{R) THEORIES IN PALATINI FORMALISM AND AUTONOMOUS EQUATIONS 



We consider the classes of theories given by the generalised action 

1 



S 



d 4 ^ 



-9 



2 k 



f{R) + C m + £, 



(1) 



where / is a differentiable function of the Ricci scalar R, C m and C r are the Lagrangians of the pressureless dust and 
radiation respectively, k — 8ttG and G is the gravitational constant. Motivated by recent observations we shall study 
the cosmological dynamics in these theories for a fiat Friedmann-Lemaitre-Robertson- Walker (FLRW) background 
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where a(t) is the scale factor and t is the cosmic time. 

As was mentioned above, in Palatini formalism the metric and the affine connections are treated as independent 
variables with respect to which the action is varied. Varying the action (JTJ) with respect to the metric g^ v gives (see 
e.g. Ref. [H) 



where F = df /OR and is given by 
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For the FLRW metric with pressureless dust and radiation we obtain the generalised Fricdmann equation 
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where a dot denotes differentiation with respect to t and H = a/ a is the Hubble parameter. Here p m and p r are the 
energy densities of matter and radiation, respectively, which satisfy 
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Contracting Eq. Q, and recalling that the trace of the radiative fluid vanishes, gives 

FR-2f = -K Pm . 

Using Eqs. ©-(H]) we obtain 
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where a prime denotes a derivative with respect to R. Combining Eqs. ((3 
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In the case of the Hilbert action with f = R, Eq. (jTUJ) reduces to the standard Friedmann equation: H 2 = n(p m +p r )/3. 

To obtain a clear picture of possible dynamical regimes admitted by these theories, we shall make a detailed study 
of a number of families of f(R) theories of the type ([T]), using phase space analysis. It is convenient to express these 
systems as autonomous systems by introducing the following dimensionless variables: 
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In terms of these variables the constraint equation (fT0|) becomes 

Kf) 



mm = i - m - ■ (13) 

Differentiating Eq. (fT0|) and using Eqs. ((5|l- (fT0|) we obtain 

2 |>=-3 + 3yi-y 2 -^-^ + ^§j. (14) 

Now using variables (|I2I) together with Eqs. (fT"3|) and (fl4|) we can derive the corresponding evolution equations for 
the variables y\ and j/ 2 thus: 



dy 
dy2 



d ^=yi[3-3y 1 +y 2 + C( J R)(l-y 1 )] , (15) 



where N = In a and 



(1V [-1 - 3yi + - C(i2) Vl ] , (16) 



r(p s = RF o (FR-2f)F'R 

~~ H(FR — f) (FR - f) (F'R - F) ' [ ' 



We also note that the following constraint equation also holds 

FR - 2/ = I -yi- V 2 
FR - f 2yi 



(18) 



which shows that R and thus C(i?) can in principle be expressed in terms of variables y\ and 2/2 ■ 

The behaviours of the variables y\ and 2/2 depend on the behaviour of the function C(R). In particular, the 
divergence of C{R) can prevent 2/1 and 2/2 from reaching equilibrium, as we shall see in the next Section. To proceed 
here we shall assume that C(R) is well behaved. In that case an important step in understanding the dynamics of 
such systems is to look at their equilibrium points/invariant sets and their stabilities. The fixed points (2/1,2/2) satisfy 
dyi/AN = = dj/2/diV. In this case (even when C{R) depends on R, but excluding the cases C(R) = —3, —4) we 
obtain the following fixed points: 

• Pr '■ (yi, 2/2) == (0, 1) , [We shall discuss this case further below.] 

• Pm- (2/1,2/2) = (0,0), 

• Pd : (2/1,2/2) = (1,0). 

If C(R) = —3, as is the case with the model f(R) oc R n (n ^ 1,2), we obtain (a) (2/1,2/2) = (0,1) and (b) 

(2/1,2/2) = (y[ C \o) where y[ c ^ is a constant. When C(R) = —4, the fixed points are given by (a) (2/1,2/2) = (0,0), (b) 

(2/1,2/2) = (1, 0) and (c) (2/1,2/2) = {y[ C \ 1 ~ Vi^)- Note that in both these cases the latter points in fact correspond to 
a line of (rather than isolated) fixed points. 

The stability of the fixed points in the 2-dimensional phase space (2/1,2/2) can then be studied by linearising the 
equations and obtaining the eigenvalues of the corresponding Jacobian matrices calculated at each equilibrium point 
[40|. Assuming that dC /dj/i and dC /dj^ remain bounded, we find the following eigenvalues for the above fixed points: 

. P r : (Ai,A 2 ) = (4 + C(i?),l), 

• P m : (Ai,A 2 ) = (3 + C(i?),-l), 

• Pd ■■ (Ai, A 2 ) = (-3 - C{R), -4 - C(R)) . 

In the following we shall also find it useful to define an effective equation of state (EOS), w c s , which is related to 
the Hubble parameter via H/H 2 = — (3/2)(l + w cS ). Using Eq. (Til} we find 

I F £ FR 

w c « = -2/1 + -2/2 + — + — -^^. (19) 

Once the fixed points of the system are obtained, one can evaluate the corresponding effective equation of state by 
using this relation. 
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III. COSMOLOGICAL DYNAMICS FOR MODELS BASED ON f(R) THEORIES 

In this section we shall use the above formalism to study a number of families of f(R) theories, recently considered 
in the literature. 

A. f(R) = R-A 

The simplest model not ruled out by observations is the ACDM model. To start with, therefore, it is instructive to 
consider this model in this context, as it represents the asymptotic state in a number of cases considered below. The 
Lagrangian in this case is given by 

f(R) = R-A, (20) 
which gives C(R) = 0. Eqs. (fT5")l and HH) then reduce to 

^ = yr(3-3y 1 + y 2 ), (21) 

^=i&(-l-3i&+itt), (22) 

with fixed points given by 

P T - Q/1,2/2) = (0,1), P m : (2/1,2/2) = (0,0), (2/1,2/2) = (1,0). (23) 

The corresponding eigenvalues can be evaluated to be 

P r - (Ai,A 2 ) = (l,4), P m : (Ai, A 2 ) = (3, —1) , P d : (Ai, A 2 ) = (-3, -4) . (24) 

Thus the fixed points P r , P m and Pd correspond to an unstable node, a saddle point and a stable node, respectively. 
From Eq. (fTT))) the effective equation of state is in this case given by w c fr = —2/1 + 2/2/3. This then leads to w c ff = 1/3, 
and —1 for the three points given in Eq. (|23p . indicating radiation-dominated, matter-dominated and de-Sitter 
phases, respectively. We have confirmed that this sequence of behaviours does indeed occur by directly solving the 
autonomous equations (f2"Tj) and (12"2")) . 

The ACDM model corresponds to f(R) = R — A, in which case C(R) vanishes because F' — 0. Equation (fTT)) 
shows that C(R) also vanishes when 

FR - 2/ = , (25) 

provided that other terms in the expression of C(R) do not exhibit divergent behaviour. When C(R) — ► the system 
possesses the fixed points P r , P m and Pd- In principle, the nature of the point Pd depends on the nature of the theory 
under study, i.e. the system (fT5|) - ([T6| . If the last terms in these equations vanish sufficiently fast as C(R) — ► 0, then 
P d corresponds to a de-Sitter solution which is a stable node since its eigenvalues are given by (Ai, X 2 ) = (—3, —4). 
If, on the other hand, these terms fall slower than 1/N, then they can contribute to the evolution of 2/1 and 2/2 and 
the point Pd may be different from a de-Sitter point. We note, however, that for all the theories considered in this 
paper Pd corresponds to a de-Sitter point. 

B. f(R) theories with the sum of power-law terms 

The families of models we shall consider in this section belong to the classes of theories given by 

f(R) = R + aR m - (3/R n , (26) 

where m and n are real constants with the same sign and a and (3 have dimensions [mass] 2 *- 1 "" 1 -* and [mass] 2 ^ n+1 ^ 
respectively. Such theories have been considered with the hope of explaining both the early and the late accelerating 
phases in the universe [23|, [HJ ■ In this subsection we shall make a detailed study of models based on such theories in 
order to determine whether they admit viable cosmological models, with both early and late time acceleration phases. 
In this case equation (p~8|) becomes 

1 - (to - 2)aR m - 1 - (n + 2) j 9fl-"~ 1 = 1 - Vl - y 2 
(m-^aR" 1 - 1 + (n + l)f3R- n - 1 2y x [ ' 
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which in principle allows R to be expressed in terms of y\ and y 2 , at least implicitly. The function C(R) is given by 

_ [1- (m-2)ai? m - 1 - (n + 2)f3R- n - 1 ][m(m-l)aR m - n(n + l)/3R- n ] 
^ ' ~ _ [1 - m(m - 2)aR m ~ 1 + n(n + 2)(3R- n - 1 ][(m - l)aR m + {n + l)/3R- n ] ' ^ ' 

Before proceeding, some comments are in order concerning the variables y\ and y 2 . From the expression for Hubble 
function (fTU|) , one can observe that if FR — f > then one requires 6F£ > 0, as otherwise H 2 would be negative. 
Assuming (see below) that at late times Eq. J2H]) tends to R — (3/R n , then FR — f — * (1 + n)[3R~ n which is thus 
positive since observations (see section TlV Ap require n > — 1 and [3 > 0. Hence at late times y\ and y 2 are positive 
and, from Eq. (|13p . their sum must be smaller than 1, as otherwise p m would be negative. Similarly assuming (see 
below) that at early times Eq. (f2"o| tends to R + aR m , then we have FR — / — > (m — \)aR m . To obtain an early time 
inflation, we need m > 1 (as we shall see below). Moreover, for the action to remain positive at early times, we require 
a > which implies / > 0. Therefore FR — f is again positive. Thus y\ and y 2 are positive also at early times, with 
a sum which is smaller than 1 for physical reasons. This indicates that at both early and late times, variables y\ and 
yi can be treated as normalised. 



Limiting behaviours 



To study the full behaviour of the theories of the type (p6|) it is useful to first consider the limiting cases where one 
of the nonlinear terms in the action will dominate. We shall consider these cases separately. 

(i) The a = case 

In this case reduces to 



and C(R) becomes 



f(R) = R — P/R n , (29) 



^)=3n£ +n -(f + ^, (30) 



which allows R to be expressed in terms of y\ and j/2 thus: 



R i+n = P [ 3 2/i + n (Vi - V2 + 1) - V2 + 1] , 31 ^ 
The critical points in this case are given by 

• Pr- (yi,y 2 ) = (0,1). 

Since the numerator and the denominator of Eq. pip tend to zero as one approaches P r , care must be taken in 
this case. We therefore split the analysis into three parts: 

1. P r i'. This point corresponds to (3/R 1+n ^ 1 and C — > in. The eigenvalues in this case are given by 
(Ai, A 2 ) = (3n + 4, 1) with w cff = 1/3. 

2. P r 2'- This point corresponds to (3/R 1+n ^> 1 and C — > — 3. The eigenvalues in this case are given by 
(Ai,A 2 ) = (1,1) with w cS = -2/3- 1/n. 

3. P r 3 : — > constant. This case, however, does not occur since then j/i ^> y 2 — 1 (the reasoning would 
still remain the same if y\ is of the same order of magnitude as 1 — J/2) which would imply that R n+1 should 
tend to /3(3 + n)/2. But from Eq. ([131), one can deduce that np m ~ (-Fi? + /)/2 and from Eq. (JSJ) that 
«/),„ = — FR + 2/. Hence this would imply (—FR + /)/2 ~ — FR + 2/ which is not possible. 

• Pm : (ift.ltt) = (0,0). 

In this case C — * 3n, (Ax, A 2 ) = (3(n + 1), —1) and w e g = 0. 

• Pd : (2/1,2/2) = (1,0). 

In this case R is a nonzero constant satisfying R 1+n = (2 + n)/3 with C = 0, (A 1? A 2 ) = (—3, —4) and w c g = — 1. 
This de-Sitter point exists for n > —2 provided that R > and /3 > 0. 
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n 


Prl 


Pr2 


Pm 


Pd 


P 


n < -2 


saddle 


unstable 


stable 


- 


saddle 


-2 < n < -4/3 


saddle 


unstable 


stable 


stable 


saddle 


-4/3 <n < -1 


unstable 


unstable 


stable 


stable 


saddle 


-1 < n < 


unstable 


unstable 


saddle 


stable 


stable 


n > 


unstable 


unstable 


saddle 


stable 


saddle 



Table I: Fixed points and their natures for models based on theories of the type f(R) = R — f3/R n . The existence of the Pd 
point here assumes that > 0. For negative the symbol '-' in the Table should be replaced by 'stable' and vice versa. 



•P- (m,y 2 ) = (-^ ) o). 

In this case C = 3n, (Ai, A 2 ) = (1 + 1/n, —1) and w e s = — 1 — l/n. 

This fixed point is, however, not relevant in the asymptotic regimes that we shall consider. This is because with 
observationally motivated parameters (see below) ^±i > o which thus implies y\ < 0. However we have shown 
that for late and early times < y\ < 1. 

A summary of fixed points in this case together with their stability properties is given in Table fl] An inspection of 
this table shows that the sequence of radiation (P r i), matter (P m ) and de-Sitter acceleration (Pd) can be realized for 
n > -1. 

(ii) The (3 = case 

In this case (121))) reduces to 

f(R) = R + aR m . (32) 

The fixed points in this model can be obtained from the previous case by simply taking n — > —m and (3 — > —a. Hence, 
depending on the energy region the fixed points can be summarised as follows. 

1. High energy points (ai? m_1 ^> 1) 

(la) P r2 : (yi,y 2 ) = (0, 1), w cS = -2/3 + 1/m and (Ai, A 2 ) = (1, 1). 

(lb) P: ( yi ,y 2 ) = (-gf^,0), w cB = -1 + 1/m and (Ai,A 2 ) = (1 - l/m,-l). 

2. Low energy points (ai? m_1 <§; 1) 

(2a) P rl ; {yi,V2) = (0,1), w cS = 1/3 and (Ai,A 2 ) = (4-3m,l). 
(2b) P m : ( yi ,y 2 ) = (0, 0), w cS = and and (A x , A 2 ) = (3(1 - m), -1). 

3. de-Sitter point (aR™- 1 = l/(m - 2)) 

(3) P d : (yi,y 2 ) = (1,0), ^ e fF = -1 and (Ai,A 2 ) = (-3,-4). 

Note that the case m = 2 requires a separate analysis since this corresponds to C = —6 independent of R [see 
Eq. (|30p]. We shall return to this case below. The fixed points and their properties for this case are summarised in 
Table HI 

One can in principle consider the fixed points in Tables U and [H] as respectively being relevant at the future 
and the past respectively. Using the above asymptotic information we shall now consider families of theories of this 
type separately. 

2. Theories of type f(R) =R- f3/R n 

In the metric approach, models based on theories of this type were considered in Refs. [HlGjI as ways of producing a 
late-time acceleration of the universe. In addition to the difficulties such theories face (as discussed in the Introduction), 
they have recently been shown to be unable to produce a standard matter-dominated era followed by an accelerated 
expanding phase for n > [25j . Using the Palatini approach, however, such theories can be shown to be able to give 
rise to such a standard matter-dominated phase, as we shall see below. 
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m 


Prl 


Pr2 


Pm 


Pd 


P 


2 < m 


saddle 


unstable 


stable 


stable 


saddle 


4/3 < m < 2 


saddle 


unstable 


stable 




saddle 


Km < 4/3 


unstable 


unstable 


stable 




saddle 


< m < 1 


unstable 


unstable 


saddle 




stable 


m < 


unstable 


unstable 


saddle 




saddle 



Table II: Fixed points and their natures for models based on theories of the type f(R) = R + aR m . The existence of the Pd 
point here assumes that a > 0. For negative a the symbol '-' in the Table should be replaced by 'stable' and vice versa. 
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Figure 1: Evolutions of the variables j/i and j/2 for the model f(R) = R — f3/R n with n = 0.02, together with the effective 
equation of state w e s. Initial conditions were chosen to be j/i = 10 -40 and 2/2 = 1 — 10~ 5 . This demonstrates that the sequence 
of (Pri) radiation-dominated (w c a = 1/3), (Pm) matter-dominated (w e g = 0) and (Pd) de-Sitter acceleration (w e s = — 1) eras 
occur for these models. 

To proceed we shall, for the sake of compatibility with the observational constraints obtained in the following 
sections, confine ourselves to the cases with n > — 1 and (3 > 0. fn such cases, Table Q] shows that P r \ is an unstable 
node (source), P m is a saddle and Pd is a stable node. Numerical simulations confirm that such theories indeed admit 
the 3 post-inflationary phases, namely radiation-dominated, matter-dominated and de-Sitter phases (see Fig.[l]). The 
point P can never be reached because, for the considered range of n, one can show that the positive variable y\ is an 
increasing function of time, i.e. dyi/diV > (45j . The point P, however, occurs for yi — and hence can never be 
reached except by taking yi < which is unphysical since it implies a negative H 2 . 

The initial value of the ratio r — yi(0)/?/2(0) plays an important role in determining the duration of the matter- 
dominated phase. This ratio is related to cosmological parameters such as the present value of VL m that will be 
constrained by observations in the next section. For the existence of a prolonged matter-dominated epoch we require 
the condition r < 1. In this case the smaller the value of r the longer will be the matter-dominated phase. For 
large enough values of r, the solutions directly approaches the stable de-Sitter point Pd after the end of the radiation- 
dominated era. 



3. Theories of type f(R) = R + aR m - f3/R n 

In this section, we shall employ the above asymptotic analyses to study the dynamics of the generalised theories of 
this type. In particular wc shall give a detailed argument which indicates that, under some general assumptions, an 
initial inflationary phase cannot be followed by a radiation phase. 

We shall proceed by looking at the early and late dynamics separately. In this connection it is important to decide 
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which phases can be considered to be produced by the early and late time approximations f(R) = R + aR m and 
f(R) = R — j3/R n respectively. We note that in the discussion below we shall at times divide the 'early' phase 
into high energy (aR" 1 ^ 1 1) and low energy (aR" 1 ^ 1 <C 1) regimes. To proceed we start by assuming that the 
early time approximation covers inflation and radiation-dominated phases and the late time approximation covers the 
matter-dominated and dark energy phases. Thus our analysis does not cover the epochs where the two nonlinear terms 
are comparable. Such epochs can occur between early and late times or at late times. Nevertheless our assumption 
regarding early times should still be valid since R should be large in such regimes. Also observations still require that 
n > — 1 and /? > 0. 

Table |TT] shows that to have an inflationary stage followed by a radiation-dominated phase we must have m > 4/3, 
since this is the only way to obtain a saddle radiation phase P r \ (which would be required to obtain an exit from this 
phase). The de-Sitter point Pd exists for m > 2, but one can not use this fixed point for inflation followed by the 
radiation era because Pd is a stable attractor which does not admit the exit from inflation. The inflationary epoch 
can only come from the point P r 2 with m > 3 or point P with m > 3/2, in order to ensure that w c g < —1/3. 

Let us consider the case m > 2. Taking into account the fact that P r 2 is in a high energy region (aR" 1 ^ 1 ^> 1), the 
possible trajectory in this case is P r 2 — » P — » Pd, since P r i and P m are in the low energy region. All points P r2 , P 
and Pd correspond to inflationary solutions for m > 3, whereas P r 2 is not an inflationary solution for 2 < to < 3. In 
neither case do we have a successful exit from the de-Sitter epoch Pd to the radiation phase. 

When 4/3 < m < 2 it is possible to have the sequence: (i) P r 2 — > P — > P m , or (ii) P r 2 — > P r i — > P m - Note that P r 2 
corresponds to a non-standard evolution with — 1/6 < w e g < 1/12. When to > 3/2 the point P leads to an accelerated 
expansion, but the case (i) is not viable because of the absence of the radiation era after inflation. Similarly, the case 
(ii) can be ruled out since it does not possess an inflationary phase. Now since both P and P r \ are saddle points, we 
may ask whether the sequence P r 2 — > P — > P r i — * Pm can occur. To study this possibility, we proceeded numerically 
and solved the autonomous equations for a range of initial conditions close to yx = and 2/2 = 1 and a range of model 
parameters. Despite extensive numerical searches we were unable to find a radiation-dominated phase between the 
de-Sitter and matter epochs. Fig. [2] gives a typical example of such simulations showing the evolution of the variables 
2/i and 2/2 together with w e R. Clearly the radiation epoch is absent between the de-Sitter and matter epochs. The 
sequence (i) P r 2 — > P — > P m is, however, realized as expected, which is finally followed by a de-Sitter universe as 
the (3/R n term becomes important. Thus the models with 3/2 < m < 2 seem to be unable to produce an initial 
inflationary phase leading to a radiation-dominated epoch. 

The above discussion indicates that the assumption to > 3/2, which is necessary for an early-time acceleration phase 
to occur, is incompatible with the subsequent radiation-dominated phase. 

To complete this analysis we also consider the range where 1<to<4/3. In this case there are two possible 
sequences: (i) P r 2 — > P — > P m and (ii) P r \ — > P m . In the case (i) the system starts from a non-standard high energy 
phase P r 2 (ai?" i_1 3> 1), which is followed by a non-inflationary phase P. Since neither an inflationary phase nor a 
radiation-dominated epoch exist, the sequence (i) is therefore not cosmologically viable. In the case (ii) the system 
starts from a low-energy radiation-dominated phase (<xR m_1 <C 1) and is followed by a matter-dominated epoch. As 
long as the term (3/R n becomes important at late times, the sequence (ii) is followed by a de-Sitter point. Thus the 
case (ii) can at least account for the last three phases required in cosmology. When < to < 1 it is also possible to 
have the sequence P r \ — > P m — > P4. 

Finally, we shall consider the case where the term aR m is comparable to R around the present epoch. This case 
is suggested by the likelihood analysis considerations of the model parameters in the next section. When < to < 1 
this case is similar to the one for the model f(R) = R — (3/R n with — 1 < n < 0. If m > 1 this corresponds to 
a non-standard cosmological evolution with aR m 3> R in the past, which can be realized for a very large coupling 
a. To complement the discussion above, we shall also briefly consider this case. To start with let us consider the 
case to ^ 2. Since ai?" l_1 3> 1 prior to the dark energy epoch, the fixed points correspond to either (i) P r 2 with 
an EOS w e s = —2/3 + 1/m, eigenvalues (Ai,A 2 ) = (1,1) or (ii) P with an EOS w e g = — 1 + 1/to, eigenvalues 
(Ai, A2) = (1 — 1/m, — 1). Since we are considering the case to > 1, we have (i) w e ff < 1/3 and (ii) w c g < 0. This 
shows that the radiation-dominated phase is not realized unless m is unity (i.e. Einstein gravity) 

We shall consider the case to = 2 in the next subsection. 



4- The special case m = 2 

To conclude our consideration of theories of the type (|26p. we shall consider the special case of these theories with 
to = 2, which have recently attracted some attention in the literature [39l ]. 



(A) Theories of the type f(R) =R + aR 2 
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Figure 2: Evolution of yi and 1/2 together with the effective equation of state w e a for the model f(R) — R + otR m — (3/R n with 
m = 1.9, n = 1, a = 1, j9 = 10" 60 . The initial conditions are yi = 1.71 x 10" 6 and y 2 = 1 - 2.09 x 10~ 6 . The system starts 
from a non-standard radiation-type phase (corresponding to the P r 2 point with w c g = —2/3 + 1/m = —0.14) followed by an 
inflationary phase (corresponding to the P point with tu e ff = — 1 + 1/m = —0.47) and then followed by a matter-dominated 
phase (corresponding to the P m point with w c g = 0). The system finally approaches a de-Sitter point Pd because of the 
presence of the f3/R n term. 

In this case we have C = —6. The critical points, eigenvalues and the effective EOS are 

. P r : ( yi ,y 2 ) = (0,1), (Ai,A 2 ) = (-2,1). 

We note that in this case w e s is undefined, as can be seen from Eq. (|55|) in the Appendix. This, however, does 
not play an important role in the reasoning below, as we shall see. 

• P m = (0,0), (Ai, A 2 ) = (-3,-1), Woff = 0. 

This matter point exists only in the region aR <C 1 because of the relation aR = 2j/i/(l — yi — 2/2)- 

• P d = (l,0), (Ai,A 2 ) = (2,3), w eS = 0. 

This point exists only in the region aR ^S> 1. 

An important feature of this model is that the fixed point P e i in this case is an unstable node mimicking a dust 
Universe. Depending upon whether the radiation (y± <C 2/2) or the dark energy (j/2 yi) dominates initially, the 
universe starts with the EOS of matter (Pd) or radiation (P r ) and ends in a matter era (P m ) which is a stable node. 
Before reaching the final attractor P m , there can exist a time N — N s such that y\ = (1 — y 2 )/3. In particular 
this always occurs when y\ > 1/3 at early times, as would be expected in a universe dominated by dark energy. In 
this case, the EOS diverges at ^V s as we see from Eq. (1531 in the Appendix. Note that when y\ = (1 — j/2 ) / 3 , the 
scalar curvature is regular (aR = 1) with finite y\ and y 2 , indicating that there is no physical singularity, despite 
the singular behaviour of the effective EOS. This can be understood by recalling that the w e s so defined is just a 
mathematical construction for f(R) theories. 

(B) Theories of the type f(R) =R+ aR 2 - j3/R 

Let us consider the case in which the terms aR 2 and f3/R are important for early and late times, respec- 
tively. Then as long as a > and (3 > 0, it is possible to have the sequence Pd — * P r — ► P m followed by a de-Sitter 
point Pd which appears because of the presence of the term (3/R. However since P r corresponds to a divergent EOS, 
we do not have a standard radiation-dominated epoch in this model. 

Let us next consider the case in which the signs of a and (3 are different, say a > and (3 < 0. We will encounter this 
situation when we carry out a likelihood analysis of model parameters in the next section. Then it can happens that 
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C(R) diverges as R approaches either (— 2/3/a) 1 / 3 or R — > y 7- 3/3. In this case the effective EOS diverges and dyi/diV 
approaches — oo. For the model parameters constrained by observations the divergence occurs at R = (— 2/3/a) 1 / 3 
before reaching i? = ^/— 3/3, since in that case (— 2/3/a) 1 / 3 > y/—3p. We also note that the radiation-dominated 
phase is absent in this case as well if we go back to the past. 

This provides an example of what can happen when C(R) diverges, which may also occur for other values of (m, n). 



C. Theories of type f(R) = R + ainR- (3 

Finally we shall consider theories of the type 

f(R)=R+ai\nR- p. (33) 

The cosmological evolution of the models based on these theories has been studied [42|, [H| and it has been claimed 
that such theories have a well-defined Newtonian limit [13] ■ For these theories Eqs. (TTTj) and (fT5|) give 

C(R) = 3a r R -« + 2al » R - 2 \, (34) 
v ; {a - a\nR + P){R + 2a) ' V ; 

and 

R- a + 2a\nR-2P 1 - y x - y 2 



a — cdn R + p 2yi 



(35) 



These imply C — > in the limit R — > oo and C — > —3 in the limit R —> 0. We also note that assuming the denominator 
in (p£4"|) is non-zero, then C = for a constant i? satisfying i? — a + 2a\n R — 2/3 = 0. Using the constraint equation 
(|3"5|) . the critical points can be found to be 



0. The eigenvalues are given by (A 1; A2) = (4, 1) with w e g = 1/3. 
—3. The eigenvalues are given by (Ai, A2) = (1, 1) with w s — > —00. 
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•0, c = - 
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(3, —1) and w e g = 0. 



From the above discussion it is clear that one can obtain the sequence of radiation-dominated (P r i, matter- 
dominated (P m )i and de-Sitter (Pd) eras. We have checked numerically that this sequence indeed occurs. Note 
that the points P r2 and P are irrelevant to realistic cosmology, since the system approaches the stable de-Sitter point 
Pd with a constant R before reaching R — > 0. 



IV. CONFRONTING f(R) GRAVITY MODELS WITH OBSERVATIONAL DATA 

The detailed phase space analysis given in the previous Section provides a clear picture of the possible dynamical 
modes of behaviour that can occur for models based on these theories. These analyses demonstrate that none of the 
theories considered here can produce all four phases required for cosmological evolution. They, however, show that a 
subset of these theories are capable of producing the last three phases, i.e. radiation-dominated, matter-dominated 
and late time accelerating phases. As was mentioned above the occurrence of these phases is a necessary but not a 
sufficient condition for cosmological viability of such theories. To be cosmologically viable, it is also necessary that 
the subset of parameters in these theories that allow these phases to exist are also compatible with observations. 
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In this Section we shall study this observational compatibility by confronting models based on these theories with 
observational data. To proceed, we rewrite equations (jHJ), ([5]) and (fTU|) . using the redshift parameter z — a /a— 1 and 
employing the expressions p m = pmail + z) 3 and p r — p r o(l + z) 4 to obtain 

FR-2J = -3H 2 Q ma {l + z) 3 , (36) 
dR 9H 2 Q m0 {l + z) 2 



dz F'R - F 

H 2 _ 3fi m o(l + z) 3 + 6fi r0 (l + zf + //iT 2 



2 



o 6 _p 



1 _,_ 9ggF'a m0 (i+ z ) 3 



2F(F'R~F) 



2 



(37) 
(38) 



where f2 m o = K/O m o/3i/ 2 and f2 r o = Kp r o/3ifg and a subscript denotes evaluation at the present time. From Eq. (|36[) 
or (|3"T|) one can express R in terms of z for a given f(R) model. Equation ([55)) can then be used to obtain the Hubble 
parameter in terms of z. 

We note that for a give n f (R) theory sourced by cold dark matter, Q m o is uniquely determined once units are 
chosen such that H = 1 [38j. If radiation is also present, then both H and f2 r o need to be specified in order to 
determine f2 m o uniquely. In the following we choose J7 r o = 5 x 10~ 5 and H = 1 (see below). 

To constrain f(R) models, we shall use three sets of data: 

• (I) The first year data set from the Supernova Legacy Survey (SNLS) 0] with 71 new supernovae below z = 1.01, 
together with another 44 low z supernovae already available, i.e. a total of 115 SNe. 

• (II) The Baryon acoustic oscillation peak (BAO) recently detected in the correlation function of luminous red 
galaxies (LRG) in the Sloan Digital Sky Survey [l(J. This peak corresponds to the first acoustic peak at 
recombination and is determined by the sound horizon. The observed scale of the peak effectively constrains 
the quantity 



^0.35 = D v (0.35)^^l = 0.469 ± 0.017 , (39) 
0.35 

where z = 0.35 is the typical LRG redshift and Dy(z) is the comoving angular diameter distance defined as 



Dv(z) = 



D M {zf 



H(z) 



1/3 



dz 

with D M = I — ■ (40) 



(III) The CMB shift parameter. This constraint is defined by 



1089 d* 



#1089 = \/n m0 H 2 J = 1.716 ± 0.062 , (41) 

which is a measure of the distance between z = and z = 1089. It relates the angular diameter distance to the 
last scattering surface with the angular scale of the first acoustic peak. An important feature of this parameter 
is that it is model independent and insensitive to perturbations |4lj . To be able to use it one has to have a 
standard matter-dominated era at decoupling. 

For the goodness of fit we employ the standard % 2 minimization, defined by 



x 2 



<A (m° bs - mf 1 ) 2 



i=i 



where n is the number of data points and m° hs and mf 1 are respectively the observed and the theoretical magnitudes 
calculated from the model. We shall marginalise the Hubble constant by defining as usual a new x 2 : 



/+oo 
e -* 2/2 dM 
-oo 



where M = M — 51n(Ho) + 25 and M is the magnitude zero point offset. After some algebraic manipulations and 
defining 



A = E K^«! i B = gmf*-Jln(A) ; a = g^ ) (42) 

p—1 1 p — 1 1 p—1 1 
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Figure 3: The 68.3% and 95.4% confidence levels for the model based on the theory f(R) = R — /3/R" , constrained by SNLS 
data only. 

we obtain 

X 2 = A-^+ln^. (43) 

Finally we note that the BAO and the CMB shift parameter do not depend on the Hubble constant. Having only one 
measure point in each case, these tests can be used in the same way as priors. 

In what follows we shall proceed to place observational constraints on the models discussed in the previous section. 
Our choice of units such that Hq = 1 is consistent with marginalisation over Hq, similar to the choice made in the 
literature [38]. 

A. Theories of the type f(R) =R- (3/R n 
We shall first consider theories of the type 

f(R) = R-/3/R n . (44) 

As we already demonstrated in the previous section, one can obtain in this case a sequence of radiation-dominated, 
matter-dominated and acceleration phases provided that n > — 1. 

Recently Amarzguioui et al. [38[ have found that models of this kind are compatible with the supernova 'Gold' data 
set from Riess et al. [6j together with the BAO constraint ([31?]) and the CMB constraint (|4"T]) , subject to constraints 
on the values of the parameters (3 and n. In their combined analysis the best-fit model was found to be (3 = 3.6 and 
n = —0.09. Here we extend this work by employing the more recent first year Super-Nova Legacy Survey (SNLS) 
data as well as taking into account the presence of radiation. In Fig. [3] we show observational contour plots at the 
68% and 95% confidence levels obtained from the SNLS data only. 

In this case the best-fit value (smallest % 2 ) is found to be \ 2 = 116.55 with (3 = 12.5 and n = 0.6. This is consistent 
with the results obtained in Ref. [HI, namely f3 = 10.0 and n = 0.51 and with the results obtained in the previous 
section according to which the transition from the matter-dominated era to the de-Sitter era occurs for n > — 1. As 
n gets closer to —1, it is difficult to realize a late-time acceleration phase since the model (l4"4l approaches Einstein 
gravity. 

If we further include constraints from BAO and CMB, the likelihood parameter space is reduced significantly as 
seen from Fig. 2J with narrower ranges of the allowed values of n and (3 lying in the intervals n £ [—0.23, 0.42] and 
(3 6 [2.73, 10.6] at the 95% confidence level. This difference is mainly due to the fact that the BAO data better 
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Figure 4: The 68.3% and 95.4% confidence levels for the model based on the theory f(R) = R — f3/R" , constrained by SNLS, 
BAO and CMB data. 
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Figure 5: Evolution of the effective equation of state w e e, for the model based on the theory f(R) = R — /3/R n , with the best-fit 
value f3 — 4.63 and n — 0.027. For very large z, w e B = 1/3, which then decreases and approaches at redshift of few tens. 
The accelerated expansion occurs around z = 1 after which u> e ff becomes smaller than —1/3. In the future w e g asymptotically 
approaches — 1. 



constrains the cold dark matter parameter than the SN data alone. The best-fit values are found to be (3 = 4.63 and 
n = 0.027 with X 2 = 116.69. 

Note that the ACDM model corresponding to j3 = 4.38 and n — lies in the 68% confidence level. In Fig.0 we have 
plotted the evolution of w e g for the best-fit case. This shows that the models in agreement with the data undergo a 
transition from a radiation-dominated to a matter-dominated epoch followed by an acceleration phase which in the 
future asymptotically approaches a de-Sitter solution. This is a nice realisation of the trajectories that follow the 
sequence P r \ —* P m — > Pd described in the previous section. 
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Figure 6: The 68.3% and 95.4% confidence levels for the model f(R) = R+ aR 2 - f3/R constrained by SNLS, BAO and CMB 
data. 

B. Theories of the type f(R) =R + aR m - f3/R n 

We next consider theories of the type 

f(R) = R + aR m - f3/R n . (45) 

In the previous section we showed that an initial inflationary phase ( requiring m > 3/2) does not seem to be compatible 
with a subsequent radiation-dominated epoch. However, dropping the requirement of an inflationary phase, such 
theories can produce the last 3 phases, namely radiation-dominated, matter-dominated and late acceleration epochs, 
provided the radiation-dominated era starts in the aR rn ^ 1 <C 1 regime. If the term aR m becomes smaller than f3/R n 
prior to the decoupling epoch, observational constraints on the models based on (|4"5|) are similar to those on the 
models based on (g4]), i.e., n E [-0.23,0.42] and (3 g [ 2.73, 10 .6] at the 2a level. 

When TJi — 2 we also showed separately in Section fill B 41 that the radiation-dominated epoch does not exist. To 
study the observational constraints in such cases, we have carried out a likelihood analysis for the model f(R) = 
R + aR 2 — (3/R considered in Ref. [39(. Note that the four parameters a, (3, m and n cannot be constrained at the 
same time with the observations used here, so we have chosen special values for m and n. With these choices, we have 
found that the model is compatible with the observational data subject to the parameter constraints a € [2.30, 3.29] 
and (3 £ [—0.012, —0.004] (see Fig. [5]). However there is no radiation-dominated stage prior to the matter-dominated 
era, as is seen in Fig. [7]which gives a plot of the evolution of w e g in the best-fit case (a = 2.69 and f3 — —0.008). 

When a > and (3 < we have already shown in Section Till B 41 that a singularity appears for w c g. In fact 
Fig. [7] shows that the de-Sitter solution is not the late-time attractor since there is a singularity in the future (around 
z —0.27) where w e g becomes (positive) infinite (and the Hubble parameter H — > 0) after a short transient period 
during which it is smaller than —1/3. 

The above discussion demonstrates the necessity of the theoretical considerations in the previous section. Using the 
observational data alone, we would obtain the misleading result that the above model is consistent with observations 
despite the absence of the radiation-dominated era. This is due to the fact that the observations we have considered 
do not probe the radiation-dominated phase. 

So far we have considered the cases with m > 1. We have also checked that the sequence of radiation-dominated, 
matter-dominated and de-Sitter eras can be realized for < m < 1, although in this case the initial inflationary 
epoch is absent. One can perform a likelihood analysis by taking several values of m and n which exist in the regions 
< m < 1 and < n < 1, while allowing a and (3 to vary. Interestingly, there are values of (m,n) in these ranges 
which fit the observational data. Table IIIII shows the best-fitting model parameters for some values of the parameters 
m and n. For this range of parameters the best fitting values of a and (3 are always found to be negative and positive, 
respectively. As n is increased we find that the magnitudes of a and f3 both increase, while the former still remaining 
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Figure 7: Evolution of ui e ff for the model R + aR 2 — (3/R with the best-fit value a = 2.69 and j3 = —0.008. There is no 
radiation-dominated epoch prior to the matter-dominated era (w e s = 0). Although w e s can be smaller than —1/3 around the 
present epoch, the late-time evolution does not correspond to a de-Sitter universe. 





m = 0.1 


m = 0.5 


m = 0.8 


n = 0.1 
n = 0.3 
n = 0.5 
n = 0.8 


(116.68,-1.33,3.37,0.27) 
(116.64,-2.47,2.54,0.27) 
(116.64,-2.85,2.56,0.27) 
(116.99,-2.85,4.71,0.28) 


(116.64,-0.22,4.66,0.27) 
(116.84, -0.45,5.95,0.27) 
(116.63,-0.82,6.36,0.27) 
(116.76,-1.03,9.16,0.28) 


(116.61,-0.08,4.90,0.27) 
(116.76,-0.24,6.15,0.27) 

(118.29,-0.47, 12.4,0.28) 



Table III: The best-fitting values of the parameters (x 2 , a, (3, fi m o) f° r the models based on theories f(R) = R + aR m — (3/R 
with some particular values of (m,n). 




negative. The increase in the amplitudes of a and (3 tends to compensate for the decrease in the 1/R n term. On the 
other hand as m increases, the amplitude of a decreases (while still remaining negative) while that of (3 increases. 
These effects have the consequence of minimising the contribution of R m relative to 1/R n . We also did not find the 
divergence of w e g in such cases. Note that from the Table UTT1 the case with n = 0.8 and m = 0.8 can produce late-time 
acceleration. This provides an example of a case where the two nonlinear terms are comparable and necessary to 
produce the acceleration at late times. There would be no late-time acceleration if a single term is chosen with such 
an exponent. 



C. Theories of the type f(R) = R + a In R - (3 
Finally we consider theories of the type 

f(R) = R + alnR- P, (46) 

where a and (3 are dimensionless constants. 

The combined constraints from the SNLS, BAO and CMB data confine the allowed values of a and j3 to the ranges 
a e [—1, 1.86] and (3 £ [—0.96, 8.16], respectively (see Fig. [5]). The best-fit model corresponds to a ~ 0.11 and (3 = 4.62 
with x 2 = 116.69. These values are not too different from the ACDM model with a = and /3 = 4.38, indicating that 
the model is not significantly preferred over the ACDM model. 

Note that the case (3 = is excluded by the data. This therefore shows that the lni? term alone cannot drive 
the late-time acceleration; one still requires a cosmological constant. Thus based on the observational constraints 
considered here, only strong theoretical reasons would motivate the consideration of theories of this type. 
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Figure 8: 68.3% and 95.4% confidence level for the theory f(R) = R + alnR- (3 constrained by SNLS, BAO and CMB data. 
The models without a cosmological constant (/3 = 0) are excluded by the data. 

V. CONCLUSIONS 

We have made a detailed and systematic study of a number of families of f(R) theories in Palatini formalism, using 
phase space analysis, extensive numerical simulations as well as observational constraints from recent data. These 
theories have been recently put forward in the literature in order to explain the different dynamical phases in the 
evolution of the universe, in particular the early and/or late acceleration epochs deduced from recent observations. 
The important advantages of theories based on Palatini formalism (as opposed to their metric counterparts) include 
the facts that they are second order, have the correct Newtonian limit, pass the solar system tests and are free from 
instabilities. This provides important motivation for their consideration. 

Considering the FLRW setting, we have expressed the evolution equations as an autonomous dynamical system 
in the form of Eqs. (| 1 5[) and (fl6|) . The ACDM model corresponds to a case where the expression C(R) given by 
Eq. (|17|) vanishes. Models, based on f(R) theories, in which the variable C(R) asymptotically approaches are able 
to give rise to a de-Sitter expansion with constant R at late times. This includes the models based on theories of 
the type (a) f(R) = R - a/R n , (b) f{R) = R + aR m - f3/R n , and (c) f(R) = R+alnR- j3. We have carried out 
a detailed analysis of the cosmological evolution for such models by studying their fixed points and their stabilities 
against perturbations. 

For models based on theories of the type f(R) = R—(3/R n , we have shown that the sequence of radiation-dominated, 
matter-dominated and de-Sitter eras is in fact realised for n > — 1. This is a stark contrast to the corresponding 
theories of metric formalism for which it is difficult to realize such a sequence (25|. 

Considering models based on theories of the type f(R) = R + aR m — f3/R n , we have found that while they are 
capable of producing early as well as late accelerating phases an early inflationary epoch does not seem to be followed 
by a standard radiation-dominated era. When m > 2, the de-Sitter point is stable, which prohibits an exit from the 
inflationary phase. Inflationary solutions can also be obtained for 3/2 < m < 2 (for a > 0), but in that case they 
directly exit into a matter-dominated phase which in turn leads to a late-time accelerating phase. For < m < 4/3. 
the radiation fixed point is unstable, which makes it possible to have the sequence of radiation-dominated, matter- 
dominated and de-Sitter phases without a preceding inflationary epoch. For 4/3 < m < 3/2, inflation is not possible 
and radiation, which is a saddle point, can only be preceeded by a non-accelerated phase such as P r 2 or an unknown 
phase not satisfiyng our late and early time assumptions. 

We have also placed observational constraints on the parameters of these models, employing the recently released 
supernovae data by the Supernova Legacy Survey as well as the baryon acoustic oscillation peak in the SDSS luminous 
red galaxy sample and the CMB shift parameter. We have found that both classes of models are in agreement with the 
observations with appropriate choices of their parameters. For models based on theories of the type f(R) = R — a/R n , 
the best fit values were found to be (3 = 4.63 and n — 0.027. This is consistent with the results obtained by Amarzguioui 
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et al. [38] who used the supernova 'Gold' data rather than the SNLS without taking into account the radiation. 

For the models based on theories of the type f(R) = R + aR m — (3/R n , it is not possible to constrain all four 
parameters simultaneously. Concentrating on the special class of theories f(R) — R + aR 2 — (3/R studied in litter- 
ature, we have found that they are compatible with the data subject to parameter constraints a G [2.30, 3.29] and 
(3 € [—0.012, —0.004]. Despite this agreement we have found, using a phase space analysis, that such cases do not 
seem to produce a radiation-dominted era prior to the matter-dominated epoch, thus making them cosmologically 
unacceptable. The reason for this apparent contradiction is that the combination of the SN, SDSS and the CMB 
shift parameter considered here does not provide constraints on early phases of the universe prior to the decoupling 
epoch. We have also checked that a sequence of radiation-dominated, matter-dominated and de-Sitter phases becomes 
possible for models based on theories of the type f(R) = R + aR m — (3/R n when < to < 1, in agreement with the 
data. 

Finally we have studied the compatibility of theories of the type f(R) = R + aln R — (5 with observations. Again we 
have found that such models are compatible with observations with appropriate choices of their parameters. However, 
the logarithmic term on its own is unable to explain the late-time acceleration consistent with observations. 

In summary, we have found that using Palatini formalism it is possible to produce models based on the classes 
of theories considered here, which possess the sequence of radiation-dominated, matter-dominated and late-time 
acceleration phases consistent with observations. However, we have been unable to find the sequence of all four 
phases required for a complete explanation of the cosmic dynamics. 
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Appendix 

In this Appendix we shall give the expression of w e g for the model f(R) = R + aR m — (3/R n . We shall write the 
various terms in the expression for w e s in Eq. (|19p . 
Using the relation ©, we obtain 

R _ 1 - (to- 2)aR m - 1 - (n + 2)f3R-' n - 1 
~HR ~ ~ 1 - to(to - 2)aR m - 1 + n(n + 2)f3R- n ~ 1 ' ^ ' 

which allows the third term in Eq. (|19p to be written as 

F R to(to- l^i?" 1 " 1 -n{n+l)l3R- n - 1 



3HF HR 3[1 + maR m ~ 1 +n/?i?-"- 1 ] 

Defining the variable £, as 

3F'(FR - 2/) 
^ = 2F(F'R - F) 

3[m(m - l)ai?'™- 1 - n(n + l^i?-"- 1 ] [-1 + (m - 2)aR m - 1 + (n + 2)f3R-' n - 1 ] 
2(1 + muR™- 1 + nfiR-™- 1 ^-! + m{m - 2)aR m ~ 1 - n(n + 2)/3i?-"- 1 ] ' 

allows the fourth term in Eq. to be expressed as 

j = R_ 

3H£ 31- (HR 

m{m - ifaR™' 1 + n(n + iffiR' 11 ^ 1 (to - 2) (to - ^aR™- 1 - (n + l)(n + 2) / 3iT"- 1 
to(to - l)ai?™- 1 - n(n + l)^""" 1 1 - (to - 2)aR m - 1 - (n + 2)/3i?-"- 1 

to(to - ^aR™- 1 - n(n + l)/?^"" 1 to(to - l)(m - 2)aR m ~ 1 + n(n + l)(n + 2)/3i?"™" 1 
1 + maR m - 1 + n/Ji?-"- 1 ' 1 - to(to - 2)aR m - 1 + n(n + 2)(3R- n ~ l 



(48) 



(49) 



(50) 
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Finally the last term in Eq. (| 19[) is given by 

FR _m{m-l)aR m - 1 -n(n+l)/3R- n - 1 R 
18F£H 3 ~ 3[(m - l)ai? m ~ 1 + {n + l^i?-™" 1 ] Vl ~HR ' ^ ^ 

Note that from Eq. (|27|) i? is a function of yi and yi, which in turn implies that w c g is also a function of y\ and j/2, 
albeit of a very complicated form. 

Despite this complexity, however, the expression for w e s reduces to a fairly simple forms for specific cases. For 
example, when m — 2 and (3 — 0, we find 

1 2aR(2 + aR) 

w cS =y 1 + -y 2 + {1 + 2aR){1 „ aR) • (52) 

Now since from Eq. (|27| aR — 2y\/ {\ — y\ — yi) in this case, w e g expressed in terms of y\ and yi thus 

, 1 , 82/1(1 — y a ) 

Wcff = 2/1 + -J/2 + TT— -rj 5 r ■ (53) 

3 (1 + 3yi - y 2 (1 - 3yi - y 2 ) 
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